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Abstract. We present initial attempts to include the multi-dimensional nature of radiation transport in hydro- 
dynamical simulations of the small-scale structure that arises from the line-driven instability in hot-star winds. 
Compared to previous ID or 2D models that assume a purely radial radiation force, we seek additionally to 
treat the lateral momentum and transport of diffuse line-radiation, initially here within a 2D context. A key 
incentive is to study the damping effect of the associated diffuse line-drag on the dynamical properties of the 
flow, focusing particularly on whether this might prevent lateral break-up of shell structures at scales near the 
lateral Sobolev angle of ca. 1°. Based on 3D linear perturbation analyses that show a viscous diffusion character 
for the damping at these scales, we first explore nonlinear simulations that cast the lateral diffuse force in the 
simple, local form of a parallel viscosity. We find, however, that the resulting strong damping of lateral velocity 
fluctuations only further isolates azimuthal zones, leading again to azimuthal incoherence down to the grid scale. 

To account then for the further effect of lateral mixing of radiation associated with the radial driving, we next 
explore models in which the radial force is azimuthally smoothed over a chosen scale, and thereby show that this 
does indeed translate to a similar scale for the resulting density and velocity structure. Accounting for both the 
lateral line-drag and the lateral mixing in a more self-consistent way thus requires a multi-ray computation of the 
radiation transport. As a first attempt, we explore further a method first proposed by Owocki (1999), which uses 
a restricted 3-ray approach that combines a radial ray with two oblique rays set to have an impact parameter 
p < R, within the stellar core. From numerical simulations with various grid resolutions (and p), we find that, 
compared to equivalent 1-ray simulations, the high-resolution 3-ray models show systematically a much higher 
lateral coherence. This first success in obtaining a lateral coherence of wind structures in physically consistent 2D 
simulations of the radiative instability motivates future development of more general multi-ray methods that can 
account for transport along directions that do not intersect the stellar core. 

Key words. Hydrodynamics - line: formation - radiative transfer - stars: atmospheres - stars: early type - stars: 
mass loss 


1. Introduction 

The driving of hot-star winds by line-scattering of the 
star’s radiation is understood to be highly unstable to 
small-scale radial perturbations (Lucy & Solomon 1971; 
MacGregor, Hartmann, & Raymond 1979; Owocki & 
Rybicki 1984, 1985). This “Line-Driven-Instability” (LDI) 
has long been supposed to be a root cause of the exten¬ 
sive small-scale stochastic wind structure that has been in¬ 
ferred from various kinds of observational diagnostics (see 
reviews by, e.g., Owocki 1994; Feldmeier & Owocki 1998; 
Dessart 2004a). For example, black troughs of saturated 
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P-Cygni line profiles are understood to arise from the net 
backscattering of multiple line resonances occurring in a 
highly non-monotonic velocity held (Lucy 1982, 1984). 
Embedded wind shocks arising from such velocity varia¬ 
tions are moreover thought to give rise to the soft X-ray 
emission often observed from such hot stars (Chlebowski 
1989; Lucy & White 1980; Feldmeier et al. 1997). An even 
more direct diagnostic for such small-scale structure is the 
explicit, small-amplitude line profile variability often de¬ 
tected in high signal-to-noise (SN ~ 100 — 1000), optical 
spectra obtained by ground-based monitoring of recom¬ 
bination emission lines in both OB (Eversberg, Lepine & 
Moffat 1998 on £ Puppis) and WR stars (Robert 1992; 
Lepine & Moffat 1999). 
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But a key limitation in developing quantitative tests 
for the association between such phenomena and the LDI 
is that, owing to the computational expense of evaluating 
non-local integrals needed for calculating the line-force, 
time-dependent dynamical simulations of the resulting 
nonlinear flow structure have for many years been lim¬ 
ited to one dimension (ID) (Owocki, Castor, & Rybicki 
1987; Feldmeier 1995; Feldmeier et al. 1997; Owocki & 
Puls 1999; Runacres & Owocki 2002). (See also Owocki 
1999 and Gomez & Williams 2003 for some previous 2D 
attempts.) These ID simulations show development of ex¬ 
tensive radial variations in both density and radial veloc¬ 
ity, but are inherently incapable of determining the de¬ 
velopment of any corresponding variations in the lateral 
direction. In recent years, we have focused efforts on deriv¬ 
ing empirical constraints on the multidimensional proper¬ 
ties (particularly lateral scale) of wind structure (Dessart 
& Owocki 2002a,b), and developing approaches for mul¬ 
tidimensional simulation of the formation of such struc¬ 
ture from the nonlinear evolution of the LDI (Dessart & 
Owocki 2003, hereafter DO-03). The present effort builds 
on the “2D-H+1D-R” (denoting 2D Hydrodynamics, but 
with only ID Radiation transport) approach of DO-03, 
now incorporating various approximate attempts to ac¬ 
count for lateral transport and radiation forces associated 
with diffuse, scattered radiation, within similar 2D hydro- 
dynamical simulation models. 

A key aspect of line-driving in hot-star winds regards 
the Doppler shift associated with the wind acceleration 
and expansion. Along any given direction n, this effec¬ 
tively desaturates the line-transfer on scales of the di¬ 
rectional Sobolev length l n = i> t h/n • V(n • v), where v 
is the velocity, and nth is the ion thermal speed (typi¬ 
cally a factor few smaller than the sound speed). In a 
smooth, supersonic wind with v nth, the Sobolev length 
l n « l 0 = rvth/v is of order Vth/v <C 1 smaller than the 
characteristic wind expansion scale r. The associated de¬ 
saturation of the line-transfer then allows one to write a 
purely local expression for the line-acceleration g in terms 
of the local projected velocity gradient, n • V(n • v), av¬ 
eraged over directions to the source radiation from the 
stellar core . 1 

This “Sobolev approximation” (Sobolev 1960) pro¬ 
vides the basic framework for the standard Castor, 
Abbott, & Klein (1975, CAK) model for steady, spheri¬ 
cally symmetric, line-driven stellar winds. Indeed, within 
this generalized, vector formulation, the purely local na¬ 
ture of the line-force makes even 3D time-dependent 
simulations computationally feasible. As such, this vec¬ 
tor Sobolev approach has been applied in numerous 


1 A further key simplification of this Sobolev approach is that 
the force associated with the diffuse, scattered component of 
radiation vanishes due to fore-aft symmetries of the Sobolev 
escape probability. As discussed below, the breaking of such 
symmetries for structures near the Sobolev scale leads to the 
diffuse “line-drag” effect first described by Lucy (1984). 


multi-dimensional models aimed at studying large-scale 
variations in line-driven winds. Associated with a rota¬ 
tional modulation of photospheric properties (Cranmer 
& Owocki 1995; Dessart 2004b), they reproduced some 
of the key features of large scale variations observed in 
blueshifted absorption troughs of UV P-Cygni profiles 
in O-star (Howarth & Prinja 1989) and Be-star (Grady, 
Bjorkman & Snow 1987) winds. Wind distortion due 
to stellar rotation (Owocki, Cranmer, & Gayley 1996; 
Petrenz & Puls 2000) also provides a key explanation for 
the observed polar-enhanced mass loss of line-driven winds 
such as those of the present day 77 Car (Smith et al. 2003). 
The radiative-braking phenomenon (Gayley, Owocki, & 
Cranmer 1997) advocates the potential of radiation to 
accelerate as well as decelerate stellar winds in massive 
binaries, a duality that proves essential to explain the ge¬ 
ometry of wind-wind collisions in massive binary systems 
(van der Hucht & Williams 1995). 

Unfortunately, such a local Sobolev approach cannot 
be used to model structure arising from the LDI, since 
this occurs at scales near and below the Sobolev length 
(Owocki & Rybicki 1984, 1985). Instead the line-force 
must be computed from a non-local radiation transport 
that can be cast approximately in terms of integral escape 
probabilities (Owocki & Puls 1996). In ID simulations, 
these escape probabilities can be evaluated from spatial 
integrals carried out along a restricted set of near radial 
directions intersecting the stellar core, repeated for a set 
of n x ~ 3voo/vth ~ 1000 frequencies that resolve the line 
thermal width over the full range of velocity shifts within 
the wind. This makes even ID simulations of the LDI quite 
computationally expensive. 

In 2D or 3D, a proper treatment of the lateral trans¬ 
port requires such integration along a more complete set of 
oblique rays ranging from transverse to radial in direction. 
A severe complication is then that nonradial integrations 
from any given grid node do not generally intersect any 
other grid nodes. As such a straightforward “long charac¬ 
teristic” approach would require repeating the integration 
anew for each grid node, with a complex interpolation for 
the variation of the flow variables along the ray. For even a 
2D grid of n r radial and azimuthal zones, this requires 
UrHcf, integrals involving of order n r n x operations for each 
of the set of n ray directions, giving an overall scaling of 
n r ayn x n x nct, operations. For a typical case with n ra y ~ 10, 

« 100 , and n x ~ n r « 1000 , this implies of order 10 12 
operations to evaluate the radiative force at each time step 
of a simulation model! 

Such timing might be reduced somewhat by a “short 
characteristic” approach that builds up the local escape 
probabilities based on the evaluation in neighboring zones 
(e.g., van Noort, Hubeny, & Lanz 2002). But before at¬ 
tempting to develop such a complex and costly general 
method, we explore here some more tractable, approx¬ 
imate treatments for the lateral radiation transport and 
associated force, aimed at gaining some initial insights into 
the key dynamical effects within 2D models. 
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In particular, a key issue in such instability simulations 
is taking proper account of the “line-drag” effect of the dif¬ 
fuse, scattered radiation (Lucy 1984; Owocki & Rybicki 
1985). In ID simulations, the associated reduction in the 
net growth rate of the strong radial instability can be mod¬ 
eled via a “Smooth Source Function” (SSF; Owocki 1991, 
Owocki & Puls 1996) method. This ignores any variations 
in the scattering source function, but accounts, through 
the radial integrations for the non-local escape probabil¬ 
ity, for small-scale, fore-aft asymmetries that give rise to 
this diffuse drag (Owocki & Rybicki 1985). This nearly 
stabilizes the flow near the wind base, but in the outer 
wind there is still a strong net radial instability from driv¬ 
ing by the direct radiation from the stellar core. 

For lateral directions not intersecting the stellar core, 
the 3D linear stability analysis by Rybicki, Owocki, & 
Castor (1990, hereafter ROC) shows that this diffuse line- 
drag leads to a strong net damping of velocity varia¬ 
tions at scales near and below the lateral Sobolev length 
lo = rvth/v. In DO-03 we speculated that such damping 
might inhibit the lateral overturning of Rayleigh-Taylor 
and thin-shell instabilities that break up radial shell struc¬ 
tures, and so might lead to an overall lateral coherence 
at an associated angular scale A</>o ss lo/r = vth/v ~ 
0.01 rad ss 1°. This “Sobolev angle” is somewhat smaller 
than, but comparable to, the typical angle scale ~ 3° in¬ 
ferred for wind structure from analysis of line-profile vari¬ 
ations in emission lines from Wolf-Rayet stars (Dessart & 
Owocki 2002a,b). 

To explore this and other effects that might influence 
the lateral scale of structure, the present paper carries 
out 2D instability simulations that account approximately 
for lateral transfer effects, including the diffuse line drag. 
Our initial approach (Sect. 2) uses a simple, local parallel 
viscosity formulation for the azimuthal diffuse line-force; 
as shown in Appendix A, the 3D linear stability analysis 
by ROC implies such a viscous scaling for the line force 
arising from azimuthal velocity perturbations on a scale 
near and above the lateral Sobolev length Iq. Results from 
our nonlinear simulations show that this lateral viscosity 
can indeed strongly damp azimuthal velocity variations, 
but does not, by itself, lead to a lateral coherence above 
the azimuthal grid scale. 

Arguing then that such coherence might instead arise 
from lateral mixing of radiation associated with the ra¬ 
dial driving, we next explore (Sect. 3) models with an az¬ 
imuthal smoothing of the radial line-force, showing that 
this does lead to a comparable lateral smoothing of the 
resulting flow structure. To account more consistently for 
both effects, we finally examine fully non-local formula¬ 
tions of both the radial and azimuthal force obtained us¬ 
ing a restricted, grid-aligned, “3-ray SSF method” first 
introduced by Owocki (1999). Initial results do show an 
extended lateral coherence of instability-generated struc¬ 
ture, but inherent limitations in the ray coverage and outer 
radial grid resolution in the method leave uncertain the 
broad applicability of these 3-ray simulations. We con¬ 


clude (Sect. 5) with a brief summary of results and their 
implications for future development of generalized multi¬ 
ray methods to account for multi-dimensional transport 
in simulations of structure arising from the LDI. 


2. Radiative Viscosity Model for the Lateral, 
Diffuse Force 


As first derived in the 3D linear stability analysis by ROC, 
a key multidimensional effect of the diffuse radiation is the 
tendency to strongly damp lateral velocity perturbations. 
Appendix A shows that this can be cast in a general an¬ 
alytic form (Eq. roii that reduces to a simple viscous 
diffusion (Eq. IQll for variations on scales near or above 
the (quite small) Sobolev length, l 0 « rv t h/v ~ O.Olr. 
To now account for the dynamical effects of diffuse lateral 
transport in nonlinear simulations of wind structure, let 
us here mimic this scaling from the linear perturbation 
analysis, and so assume that the azimuthal component of 
the diffuse line-force can be approximated as a standard 
parallel viscosity term, 


f/diff,(/l>(c) — QVis(UhO 


d 2 v <t ,(r) 

r 2 d(f) 2 


( 1 ) 


Here a v i s is a dimensionless parameter to control the over¬ 
all strength of the viscous dissipation in terms of the as¬ 
sumed dimensional scale uth r for the kinematic viscosity. 
This is indeed the scaling derived from the linear pertur¬ 
bation analysis (cf. Eq. roii . with a V i S ss sa/6 a fraction 
of order a few tenths. 

More generally, this approach is also consistent with 
the notion that any moderate- to large-scale net asymme¬ 
try in the diffuse radiation field should be set by depar¬ 
tures from the Sobolev limit defined by v^/v —> 0. In par¬ 
ticular, following the standard Sobolev approximation to 
next highest order in this small parameter leads to correc¬ 
tions that scale with Vth times the gradient of the Sobolev 
optical depth r ~ KpVth/v', where v' is the velocity gradi¬ 
ent along the direction of interest. For lateral components 
in which the gradient operator V ~ d/dq !>, the combined 
derivatives lead to terms that scale as Vthd 2 v ( j > /d(j) 2 , as in 
Eq. (fH>. 

Using this formulation, we have performed simulations 
analogous to the 2D-H+1D-R models described in DO-03, 
but now including, for negligible additional CPU costs, 
this viscous approximation m for the azimuthal diffuse 
radiation force, assuming various values of the dimension¬ 
less coefficient a v is- The radial and azimuthal grids are as 
defined in DO-03, with n$ = 60 azimuthal zones of con¬ 
stant angular size 6<j> = 0.2°, thus extending over a wedge 
of A <j> = n<f,x 6(f> = 12°, with periodic boundary conditions 
in azimuth. 

As in the 2D-H+1D-R case, the simulations begin from 
a smooth, relaxed CAK model, from which there is ini¬ 
tial formation of laterally coherent shell structures that 
arise from the strong radial instability of the line driv¬ 
ing. Over time, these shells again break up from thin-shell 
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Fig. 1. Grayscale representation of the density contrast (normalized to the initial relaxed CAK model obtained with 
identical wind and stellar parameters) for a sequence of models with distinct viscosity amplitude a v is (see Sect. 1 and 
Appendix), increasing clockwise over the range 0.01, 0.1 and 1. For readability, each computed wedge is duplicated 
and stacked four times in azimuth, using the periodic boundary conditions employed in that direction. 


and Rayleigh-Taylor instabilities, but with a final asymp¬ 
totic form that is quite different depending on the assumed 
value of a v is- Figure 1 compares a representative snapshot 
of this asymptotic density structure (normalized by the 
density in the smooth CAK initial condition) for three 
models, divided into three display wedges of 60° (rep¬ 
resenting 5 repetitions of the 12° model computational 
wedges), corresponding in clockwise order to models with 
avis = 0.01, 0.1, and 1. Comparison with the top panel of 
Fig. 1 in DO-03 shows that the structure in the low vis¬ 
cosity case a v j s = 0.01 (leftmost wedge of Fig. 1 here) is 
quite similar to that obtained in this previous model with¬ 
out any lateral radiation forces; but the more viscous cases 
with a v j s = 0.1 and a v is = 1 (middle and right wedges of 
Fig. 1) are quite different, with progressively less radial 
elongation in filamentary structures. 

To provide a more quantitative comparison, Fig. 2 
shows representative statistical properties derived from 
time- and angle-averages of the density and velocity fields. 
For the case where avis = 0.01 (left wedge in Fig. 1 and 
top row in Fig. 2), the amplitude of azimuthal velocity 
variations is slightly reduced, but otherwise the clump¬ 
ing factor /ci, the radial velocity dispersion Aiy irms , and 
the velocity-density correlation function are all similar 
to those given in DO-03 with no lateral viscous term 
(o^vis = 0). 


By contrast, for the stronger viscosity cases a v is = 0.1 
and avis = 1 (middle and lower panels of Fig. 2), we see 
that the azimuthal velocity is markedly reduced, with as¬ 
sociated changes also in the other statistical parameters. 
The inhibition of lateral motion in effect tends to isolate 
further each of the azimuthal coordinates. In this situa¬ 
tion, when a shock occurs, material from below ramming 
into the dense structure is prevented from circumventing it 
by the stronger radiative viscosity in the lateral direction. 
This behavior becomes even more pronounced in the high¬ 
est viscosity case (a v is = 1), for which the dense structures 
resemble both radially and azimuthally confined clumps, 
i.e. dots in 2D. As shown in the bottom row of Fig. 2, 
the azimuthal velocity dispersion is then reduced by a 
factor of 20-30 compared to the first case. This corre¬ 
sponds to about a tenth of the sound speed, i.e. three 
orders of magnitude less than the wind flow speed. We 
see also that the clumping and radial velocity dispersion 
increase significantly with a v ; s , converging for high a v i s 
to values similar to those found in purely ID non-Sobolev 
simulations. Then, the lateral communication becomes so 
inhibited that each direction is essentially sheltered from 
its neighbors, and the 2D computed grid merely looks like 
a series of ID non-Sobolev simulations stacked together in 
azimuth. 

For a characteristic wind radius r = 61?*, Fig. 3 com¬ 
pares the lateral variation of the azimuthal velocity and 
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Fig. 2. Top row Radial variation of time- and angle-averages that characterize the nature of flow structure for the 
case with a v is = 0.01. a) clumping factor / c i, b) radial and c) azimuthal velocity dispersion, and d) velocity-density 
correlation coefficient (Runacres & Owocki 2002). Middle row Same as top for a model with a v i s = 0.1. Bottom 
row Same as top for a model with a v ; s = 1. For ease of comparison, corresponding quantities are plotted over the 
same ordinate range. 


lateral radiation force for these three values of a v is- Note 
that the magnitude of the viscous forces are comparable 
(solid curves), even though the associated azimuthal veloc¬ 
ity (dotted curves) is much smaller for the higher viscous 
coefficient. Note also that the sign of the viscous force 
depends on the concavity of the velocity variation. 

A key result here is that the lateral diffuse radiation 
has a net effect that, in some sense, is opposite of what was 
anticipated, e.g., in DO-03. In particular, because the line- 
drag of diffuse radiation strongly damps lateral velocity 
variations on scales smaller than the lateral Sobolev length 
l Q « rvth/v, DO-03 speculated that this could inhibit the 
operation of Rayleigh-Taylor or thin-shell instabilities on 
this scale, and thus lead to a finite azimuthal coherence at 
the associated angle scale A^sob ~lo/r = v t h/v. However, 
the above results show that such diffuse radiation damp¬ 
ing, as modeled here in terms of a lateral viscosity, tends 
instead to limit further the lateral communication between 


neighboring azimuthal zones, and thus lead to an even 
greater level of lateral incoherence. 

Since lateral line-drag does not limit structure to a fi¬ 
nite azimuthal scale, we thus next consider how this might 
instead arise from the angle averaging of the backscattered 
radiation, which we next model in terms of an azimuthal 
averaging of the radial driving force. 

3. Azimuthal Averaging of the Radial Radiation 

Force 

A key limitation of the above lateral viscosity approach is 
that it still does not account for ways in which the lateral 
radiation transport might alter the radial component of 
the line-driving force. Indeed, since the radial driving at 
each azimuth is solved independently by strictly radial in¬ 
tegrations of the line optical depth, there arises, through 
the diffuse component of the radial force, a backscattering 
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Fig. 3. Left Model with a v ; s = 0.01. Azimuthal variation of the normalized viscous diffuse force gdiS^ir, <f>)/ g*(r), 
at r = 61?*, and corresponding to the time-snapshot of Fig. 1. The lateral velocity, normalized to the thermal speed, 
is over-plotted (dotted line). Middle Same as left for a model with a v j s = 0.1. Right Same as left for a model with 

Ovis = 1 * 


feedback between the outer and inner wind that is effec¬ 
tively isolated and independent for each azimuthal zone. 
As such, both the previous 2D-H+1D-R and the present 
lateral viscosity models tend to develop structure that is 
independent for each azimuthal zone, implying a lateral 
incoherence down to the azimuthal grid scale. 

More realistically, such radial backscattering should in¬ 
corporate a nonzero level of lateral averaging, associated 
with transport along oblique rays that couple zones of 
different radius and azimuth. Unfortunately, as discussed 
in the Introduction, there are severe computational chal¬ 
lenges to accounting for such oblique transport through 
direct integration along non-radial rays. 

Thus, as a first approximate attempt to explore such 
effects, we have carried out here simulations in which 
the radial force is simply averaged over azimuth using a 
Gaussian filter with a tunable “smoothing angle” a At 
a given azimuthal coordinate </>j and radial coordinate r, 
the smoothed radial force takes the form, 


Hr, smooth (U ) 


Sfcgr(r, (t> k )e-^-^/^ 2 


Figure 4 compares the normalized density structure for 
simulations with tr^ of 0, 0.1, 0.5, 1 and 1.5°, represented 
respectively by the clockwise sequence of wedges of width 
of 36°, each of which corresponding to 3 azimuthal periods 
of A<f> = 12° for each of the 5 model cases. The upper pan¬ 
els show results for the full simulated grid (out to 101?*), 
while the lower panels zoom in on the region between the 
photosphere and 31?*. Note that cases with a smoothing 
angle less than the lateral grid resolution 5(f> = 0.2° 
show little departure from standard 2D-H+1D-R simula¬ 
tions. However, as increases above 5(j>, the coherence of 
wind structures becomes more pronounced. Indeed, for a 
smoothing angle of = 3°, or twice the maximum value 
shown in Fig. 4, we have found the lateral coherence ex¬ 
tends across the full azimuthal wedge of 12°, so that the 


model effectively recovers the radially confined concentric 
shells corresponding to purely ID (spherically-symmetric) 
instability simulations. Overall, this experimentation, al¬ 
though artificial, highlights the potential role of azimuthal 
averaging of radial driving in setting the lateral coherence 
of radiatively driven wind structures. 

4. Simulations Using a 3-Ray Method for the 

Nonlocal 2D Line-Force 

4.1. Method Formulation 

While the above models use a quite intricate Smooth- 
Source-Function (SSF) escape-integral formulation for the 
radial force along each azimuth, their account of lateral 
transport effects on the radiative force is only phenomeno¬ 
logical, through either a local viscosity or a simple az¬ 
imuthal smoothing. A more consistent approach would 
carry out similar escape integrals along an appropriate 
set of oblique rays spanning the range from radial to lat¬ 
eral directions. Unfortunately, as already noted in Sect. 
1 , such a calculation presents severe computational chal¬ 
lenges, stemming largely from the general misalignment of 
these rays with the nodes of the computational grid. 

However, as first introduced by Owocki (1999), there 
is a specialized spatial grid that can allow a tractable, 3- 
ray, non-local formulation of both the radial and lateral 
components of the line-force. At any given grid point, the 
non-local escape probabilities are evaluated along one ra¬ 
dial ray, plus two nonradial rays on opposite sides of the 
radial direction, set always to have a fixed impact parame¬ 
ter p < i?* toward the stellar core. A key trick is to choose 
the radial spacing so that each nonradial ray intersecting 
a grid point with indices {i,j} will also intersect other 
grid points {i ± n,j ± n}, for integer n > 1. This avoids 
the need to carry out a conceptually complex and com¬ 
putationally costly interpolation between a ( p , z) ray grid 
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Fig. 4. Grayscale images of the density contrast normalized to an relaxed CAK model with identical wind and stellar 
parameters. At the top (bottom), models are shown only out to 10 (3) R*. Each panel is composed of five different 
models (duplicated and stacked two times laterally to improve the visibility), ordered clockwise with increasing a</,, 
i.e. 0, 0.1, 0.5, 1 and 1.5°. 


for the radiation transport, and the (r, <p) grid for the hy¬ 
drodynamics. For uniform azimuthal spacing <5</>, such ray 
alignment occurs for radial grids satisfying 


1 cos [iScj) + arccos(p/.R*)] 

Figure 5 illustrates the grid and ray alignment for 5(j> = 
0.9°, p/R* = \/079 and n r = 74. Through spatial inte¬ 
gration along the 3 such rays for each of the azimuthal 
zones, one obtains a “6-stream” description (i.e. in 2 direc¬ 
tions along each of the 3 rays) for the required non-local 
escape probabilities from each of the n r x n$ grid nodes. 

Because each ray integration provides the escape prob¬ 
abilities along ca. n r grid nodes, the overall timing is now 
linear (not quadratic, see Sect. 1) in n r , scaling overall as 


3 n<t,n r n x s=s 10 s operations, or ca. a factor 10 4 savings over 
the ca. 10 12 operations estimated for a “brute-force”, long 
characteristic computation in a non-aligned grid (cf. Sect. 

I)- 

The difference in escape along the two nonradial rays 
provides a rough treatment of the lateral radiation trans¬ 
port. Because these rays are restricted to always impact 
within the stellar core radius (p < i?*), they are best 
suited for approximating the direct component of the line- 
force; but in the crucial wind-acceleration region near the 
star, they have a substantial azimuthal component, and so 
also provide a rough approximation of the azimuthal part 
of the diffuse line-force, including, for example, the impor¬ 
tant lateral “line-drag” effect that is predicted to strongly 
damp small-scale azimuthal velocity variations (ROC). As 
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Fig. 5. Illustration of 3-rays with p/7?» = 0, ±v0l) 
spaced azimuthally every Scj = 0.9°, together with the 
radial grid spacing (circular arcs) that aligns the rays to 
intersect multiple (r, cj) grid nodes; this allows very ef¬ 
ficient evaluation of the ray escape integrals needed for 
computation of the SSF line-force in 2D wind models. For 
clarity, only every ninth azimuthal or radial zone is shown. 


the rays become increasingly radial at larger radii, this ca¬ 
pacity to approximate the lateral, diffuse radiation is lost, 
but the 3 rays still provide a quite accurate representation 
of the finite-disk form for the direct line-force. 

A more serious limitation arises from the severe loss 
of radial resolution at large radii, as demonstrated by the 
radial/azimuthal grid aspect ratio, 


Srj 

TiS(j 


V(i~i/p) 2 - 1 , 


( 4 ) 


which increases as Ti/p at large radii. This means that 
small-scale radial structure can be relatively well resolved 
in the inner wind, but then becomes increasingly more 
damped by grid averaging in the outer wind. 

In our general scenario of self-excited instability, the 
nonlinear structure at larger radii plays a key role in seed¬ 
ing perturbations in the inner wind, through the backscat- 
tering feedback of the diffuse line-force. Thus it is impor¬ 
tant to minimize as much as possible this level of outer 
grid damping, by choosing as fine as possible radial grid, 
which in the 3-ray formalism here requires using as small 
as possible azimuthal grid size 5(f>. In the simulations de¬ 
scribed below, we thus choose Scj « 10 _3 rad « 0.06° (see 
Table 1) that are even finer than in previous models. For 


Table 1. Grid parameters used in the 3-ray SSF models. 
All simulations use 60 azimuthal zones and cover the first 
10 stellar radii. We also give the extrema of the radial and 
azimuthal resolutions, to compare with the Sobolev length 
at mid wind-heights of ca. 0.01 R* along the three rays. 


Model 

6r/R * 

min. max. 

n r 

P/R * 

Scj (rad) 

A 

2.2 10" 4 

0.074 

1737 

^09 

6.65 10" 4 

B 

6.7 1CT 4 

0.101 

1076 

V^5 

6.65 10” 4 

C 

4.5 10" 4 

0.155 

869 

V^9 

1.33 10" 3 

D 

1.3 1CT 3 

0.220 

539 

V^5 

1.33 10" 3 


a characteristic wind Sobolev length Iq « rvth/v ~ O.Olr, 
this implies the radial grid size Sr is sufficient to re¬ 
solve unstable structure near this scale out to a radius 
r ~ p/lOOScj « 10 p « 71?*. 

4.2. Linear Perturbation Test for 3-Ray Azimuthal 
Diffuse Force 

Before applying this 3-ray model in simulations of the non¬ 
linear evolution of wind structure, it is instructive to ex¬ 
amine its response to linear, test perturbations applied 
in a smooth, CAK initial background model. For this, 
let us compute the diffuse force for the simple case of 
a small amplitude lateral velocity perturbation v$ that 
has a Gaussian variation in both radius and azimuth, cen¬ 
tered on a grid coordinate (r c , <j> c ), with the corresponding 
widths A r c and A <f> c , 

v^,{r,<t>) = v t h exp 

Specifically, let us choose a representative radius r c = 
1.57?*, with cj c set at the mid-angle of the azimuthal grid, 
and then consider two cases distinguished by the relative 
scale of the assumed widths. For both cases, the 3-ray grid 
parameters p and Scj) are taken from the set B in Table 1. 

The first case assumes a very small scale for both 
the radial and azimuthal widths, Ar c = 0.00057?* and 
Acjo = 0.0001 rad, smaller in fact than the corresponding 
grid sizes, so that in effect the perturbation is essentially 
confined to the single grid point at (r c ,<j c ). Since both 
scales are thus also much smaller than the characteristic 
Sobolev length l 0 ss rv th/'C ~ O.Olr, the response should 
follow the direct damping scaling given by Eq. (IA.41) of 
the linear perturbation analysis in Appendix A. The lower 
panel of Fig. 6 shows just this type of scaling for this small- 
width case of the 3-ray numerical model, with the resulting 
azimuthal force (solid curve) varying with the negative of 
the azimuthal velocity perturbation (dashed curve). 

The second case assumes a very large radial extent 
Ar c = 0.27?*, and an azimuthal extent A<fi c = 0.008 rad 
that corresponds roughly with the characteristic Sobolev 
angle lo/r « 0.01. With such comparatively large scales, 
the linear perturbation analysis in Appendix A now pre- 
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diets the response to have the diffusive scaling (Eq. roi 
which implies a force that varies as the second derivative 
of the applied Gaussian velocity perturbation. The upper 
panel of Fig. 6 shows that the 3-ray response in this larger- 
width case does indeed have just this type of scaling. 

Overall, these tests demonstrate that, despite the in¬ 
herently restricted nature of the angle quadrature, the 3- 
ray method does properly reproduce the appropriate scal¬ 
ings for the azimuthal component of the diffuse line-force. 
We thus now proceed to apply this method in simulations 
of the nonlinear evolution of the 2D flow structure. 

4.3. 3-Ray Nonlinear Simulation Results 

We have used this 3-ray method to carry out 2D simu¬ 
lations for the nonlinear evolution of wind structure for 
four different combinations of the two basic grid param¬ 
eters, namely the azimuthal resolution S(f> and the im¬ 
pact parameter p of nonradial rays, with values summa¬ 
rized in Table 1. Both choices of the azimuthal grid size, 
5(j) = 1.33 x 10 -3 and 6.65 x 10 -4 rad, are fine enough 
to well resolve the characteristic lateral Sobolev angle, 
Zo/r « 0.01. The fine angle resolution is chosen in part 
because in this 3-ray method, a smaller S(f> implies, for 
any given p , a finer radial resolution Sr (Eq. • A larger 
p also implies a finer radial grid (cf. Eq. , as well as a 
more azimuthal orientation of the nonradial rays, espe¬ 
cially near the stellar surface; this gives a greater sensi¬ 
tivity to lateral variations. The initial investigations by 
Owocki (1999) found the relaxed wind properties could 
depend quite strongly on the value of p. Our two choices 
of p = VOA and p = \/(L9 allow us to further investigate 
this sensitivity. 

In addition to the basic parameters, Table 1 also sum¬ 
marizes the number of grid points needed to reach the 
fixed outer radius R max = 101?*, as well as the extrema 
of the radial resolution for each model. The latter reveal 
that the ray-projected Sobolev length is very well resolved 
at low heights, but unresolved near the outer grid radius. 
Apart from these differing grid properties, all four models 
have the same wind and stellar parameters as the models 
presented in DO-03, and in the previous two sections. To 
isolate physical effects that result solely from the use of a 
3-ray radiation force, we also run all four models using just 
the radial ray transport; these thus represent 2D-H+1D- 
R models, but run on the (quite distinctive) 3-ray spatial 
grids. 

Figure 7 shows a montage of grayscale images of the 
relaxed wind density for the four models, each normal¬ 
ized to the starting conditions obtained with an equivalent 
CAK model. The top row shows the full radial extent (to 
Rmax. = 10i?*), while the bottom row zooms into the inner 
wind region (to 3/?*). The left (right) column corresponds 
to 3-ray (1-ray) simulations. Each panel is divided into 
four wedges, corresponding respectively to models A, B, 
C and D stacked clockwise from the vertical, i.e. in order 
of decreasing radial resolution. 



0 (deg) 


Fig. 6. Projected diffuse acceleration in the lateral di¬ 
rection from the two lateral rays (solid line) for a long- 
(top) and short- (bottom) wavelength velocity perturba¬ 
tion (dotted line). For both cases, we use a relaxed CAK 
radial velocity. 

For the lower resolution models (C and D), note that 
the shell structure is never really disrupted, irrespective 
of the impact parameter p , and, moreover, irrespective of 
the inclusion or neglect of lateral transport. In contrast, 
for the higher resolution models (A and B), the role of the 
diffuse force is quite apparent from a comparison of the left 
(3-ray) and right (1-ray) columns, especially if one focuses 
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on the bottom-row panels that concentrate on the more 
well-resolved, inner wind. Most notably, in this region the 
lateral coherence is much greater when lateral transport is 
taken into account. However, there does not seem to be a 
significant dependence on the choice of impact parameter, 
with shell break-up starting at similar heights of ca. 2 i?* 
in both models A and B. 

We emphasize that, because of the degraded radial res¬ 
olution and increasingly radial orientation of the oblique 
rays, the structure forms seen in the outer regions of these 
models are probably not very realistic. However, the inner 
regions have high resolution in both angle and radius, and 
thus we believe it is significant to find that the models 
with lateral transport can retain a lateral coherence that 
extends well beyond the azimuthal grid scale. This is quite 
distinct from the results of both 2D-H+1D-R models (cf. 
DO-03 Fig. 1 and the left and right columns of Fig. 7 here) 
and the models of Sect. 2, which use a viscous form for the 
lateral diffuse force. These results thus encourage further 
investigation of the effect of lateral transport in setting a 
finite lateral scale of instability-generated flow structure 
in line-driven stellar winds. 

5. Conclusions, Summary, and Future Outlook 

This paper is a continuation of our ongoing efforts to de¬ 
velop multi-dimensional models for the nonlinear evolu¬ 
tion of structure arising from the line-driven instability in 
hot-star winds. Building upon the previous 2D-H+1D-R 
simulations of DO-03, which carry out 2D hydrodynamics 
but use a ID, radial form for the line-transfer and force, 
we develop here approximate methods to account for two 
key effects expected in multi-dimensional treatments of 
the line-force, namely lateral line-drag and lateral averag¬ 
ing. 

A summary of key results is as follows: 

1. The lateral drag effect of diffuse line-radiation can be 
well modeled by casting the azimuthal line-force in a 
parallel viscosity form that scales with the second spa¬ 
tial derivative of the azimuthal velocity. For viscous 
coefficients of order the values derived from linear per¬ 
turbation analysis, this does indeed lead to a strong 
damping of lateral velocity variations in 2D, nonlinear 
simulations of the line-driven instability. 

2. However, contrary to previous expectations, this re¬ 
duction of lateral flow does not lead to a coherence at 
a characteristic Sobolev angle (ca. 1°) associated with 
the strongest damping. Instead, it tends to further iso¬ 
late the azimuthal zones of the model, and so again 
leads to strong lateral incoherence down to the grid 
scale. 

3. Lateral coherence may instead arise from the lateral 
averaging associated with the radial driving, as demon¬ 
strated here from heuristic models that explicitly im¬ 
part an azimuthal smoothing to the radial driving 
force. 


4. To account more consistently for both the lateral line- 
drag and lateral averaging effects, we have also exper¬ 
imented with a restricted, grid-aligned, 3-ray method 
that provides an approximate nonlocal treatment for 
both the radial and lateral components of the trans¬ 
port and line-force. Results for the highest resolution 
simulations indicate lateral averaging does tend to sus¬ 
tain an extended lateral scale in the resulting flow 
structure, but specifics of the results are left uncertain 
by the inherent limitations of the method. 

Overall, a central conclusion is thus that lateral aver¬ 
aging of diffusion radiation appears promising as an effect 
that could set a distinct lateral scale to structures arising 
from the line-driven instability. However, further model¬ 
ing with a less-restricted multi-ray method will be needed 
to confirm this idea, and to determine quantitatively the 
likely scale. In particular, as mentioned in the introduc¬ 
tion and in DO-03, Rayleigh-Taylor or thin-shell insta¬ 
bility likely plays a competing role in controlling such a 
lateral scale, although with a magnitude difficult to assess 
at present. We intend to explore development and appli¬ 
cation of such methods in our continuing work within this 
overall effort, with results to be reported in future papers 
in this series. 

Acknowledgements. SPO acknowledges support of NSF grant 
AST-0097983, awarded to the University of Delaware. 


Appendix A: Linear Analysis for Lateral Diffuse 
Force 

Some key insights into the role of the multidimensional diffuse 
line force can be gleaned from the 3D linear perturbation anal¬ 
ysis carried out by ROC. Let us begin with ROC’s Eq. (A2) 
for the perturbed radiative acceleration tensor with compo¬ 
nents Tij = Sgi/5vj, giving the i’th component of the radiative 
acceleration Sgi arising from the j’th component of the velocity 
perturbation 5vj. For the diffuse component (stemming from 
the second term in the square bracket), we first note from sym¬ 
metry and parity arguments that the tensor is purely real and 
diagonal, given by (for the pure scattering case e = 0) 


Ta — sH 


2 (nik/2x*Qo) 1 2 3 \ 
n * 1 + ( mk/2x*Q 0 ) 2 / 


(A.l) 


where the angle brackets denote averaging over solid angle, and 
we have assumed the usual case of longitudinal perturbations 
with k || Sv. The factor Q o, defined in ROC Eq. (13), accounts 
for the angle variations of the velocity gradient in the smooth, 
spherically symmetric background flow, but for simplicity we 
have approximated this here by its value Qo = v 0 / (ruth) = 1 /Zo 
at the isotropic expansion point, where o = dlnv 0 /d\nr — 1 = 
0. This makes the order unity factor s = (D) / (Dfi) = 2/(1 + 
/r*), where /r* = yf% — Ri/r 2 . The overall factor — sD gives a 
net damping rate that scales with the radial instability growth 
rate, Q = 2x»ag 0 /veh, where a « 0.6 is the CAK exponent. 
For simplicity we henceforth take the blue-edge frequency x * 
to be unity. 
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Fig. 7. Grayscale representation of the density contrast for relaxed non-Sobolev models computed with the 3-ray (left 
column) and 1-ray (right column) approaches. All models were run up until a total real time of 120,000 seconds. Panels 
in the top (bottom) row show models only out to 10 (3) R*. Each panel is sub-divided into four wedges corresponding 
to models A, B C and D, stacked clockwise from the vertical in order of decreasing radial resolution (Table. 1). As 
before, each computed angular wedge is duplicated a number of times so that the total angular extent of each wedge 
shown equals 22.5°, i.e. a quarter of the quadrant in each of the four panels. The key result shown here is the presence 
of laterally coherent wind structures close to the photosphere in the high resolution 3-ray models (A and B), while 
with the equivalent 1-ray model, one observes instead wind structures with no noticeable lateral coherence, i.e. of the 
order of the lateral numerical grid resolution. 


In this work, we are interested in the lateral components ordinates (0,<j>) takes the form m = n^, — sin 6 cos <j>. Then 
Tn = T 22 , representing variations along the azimuthal direc¬ 
tion, for which the direction cosine in standard spherical co- 
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defining I\ = k/2Qo, the integrals required for evaluation of 
the angle averaging take the form, 


5g<t> 

5vd, 


sLl 

47T 


d(f) / dg 


K 2 cosV(l -M 2 ) 2 
1 + A' 2 cos 2 </(l ’ 


(A.2) 


where fi = cos 9. From numerical evaluation of these integrals, 
we find that this can be approximated (within ca. 10%) by the 
simple form, 


fig,/, ^ sQ K 2 sfi k 2 

5^ ~ 3” 1 + K 2 = 3~ 4 Ql + k 2 ' 


In the short-wavelength limit k 2Qo, this gives 


dg^ 0 _ sQ 2 asg D 

~ J ldamp — 0 — 0 

ovcf, 3 3-uth 


k > 2Q 0 , 


(A.4) 


where Lldamp represents a damping rate associated with the 
lateral line-drag of this diffuse radiation. 

In the long-wavelength limit k <C 2Qo, we find 


^ « -(as/6) (u th r) k 2 ; k « 2Q 0 , (A.5) 

4y 0 

where the latter equality uses the approximation g 0 w v^/r. 
Note that this is just the form that would arise from a parallel 
viscosity with kinematic viscosity coefficient proportional to 
Vth r, as assumed in Eq. CD- This thus forms a key motivation 
for the lateral viscosity approach used in Sect. 2 to model the 
2D nonlinear development of structure arising from the strong 
radial instability of line driving. 

The limiting forms 1A.4I and llA.511 of Eq. 1A.3I also pro¬ 
vide a means to understand the small- and large-scale linear 
perturbation tests done for the 3-ray SSF method in Sect. 
4.2. For the case in which the Gaussian perturbation width 
rA(j) 0 = 0.0005 r is much smaller than the lateral Sobolev 
length lo = rvth/vo ~ 0.01 r, the lower panel of Fig. 6 shows 
that g<f, m — 0.5g*Vrt,/vth, which for a characteristic CAK ac¬ 
celeration g 0 ~ g*/{ 1 — a) is quite consistent with the scaling 
given in Eq. <A.4t for this small-scale limit. 

However, for the larger scale perturbation, the upper panel 
of Fig. 6 shows that the force response scales as the second 
derivative of the original Gaussian velocity perturbation, con¬ 
sistent with the viscous scaling predicted in the large-scale limit 
of Eq. IIA.5I . 
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